/*!  \file  GK.h
 *	  \brief class that reporesents  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      09/02/2013 11:06:06 AM
 *   \copyright GNU Public License.
 */
#pragma once
#include "macros.h"
#include "CompileTime.h"
#define EIGEN_FUNC_PARAMS double * RS x, const int incX, double & xNorm2, double & residual, double normUV[2]

class Matrix;
class GK;
class GKManager;
typedef void (GK::*factorize)(const GK & mat, const double shift);
typedef void (GK::*initializeF)(const double shift);
typedef void (GK::*printF)(std::ostream & os ) const;
typedef int  (GK::*countF)(const double s ) const;
typedef void (GK::*addShiftF)(const double & d ) ;
typedef void (GK::*solveCanonicalVector)(const int k, EIGEN_FUNC_PARAMS) const;
typedef void (GK::*solveVecor )(const double * RS b, const int incB, double * RS x, const int incX) const;
typedef void (GK::*getEigenVec)(EIGEN_FUNC_PARAMS) const; 
typedef double (GK::*checkRepresentationF)(int &) const; 

struct Interval {
	double center;
	double width;
	int  leftCount;
	int rightCount;
	int count() const { return rightCount - leftCount; };
	void getLeftRight(double & left, double & right) const {
		left  = center - width/2;
		right = center + width/2;
	}
	void update(const int countCenter, double & left, double & right) {
		left  = countCenter == leftCount  ? center :  left;
		right = countCenter == rightCount ? center : right;
		center = (left + right)/2;
		width  = right - left;
	}
	void split(const int countCenter, double & left, double & right, Interval & oth ) {
		oth.center     = (right + center) / 2;
		oth.width      =  right - center     ;
		oth.leftCount  = countCenter;
		oth.rightCount = rightCount;
		right 			= center;
		center     		= (right + left) / 2;
		width      		=  right - left ;
		rightCount 		= countCenter;
	}
	bool isClustered(const double gapTol, const Interval & next) const {
		assert(rightCount == next.leftCount);
		double right = (next.center - next.width/2);
		double left  = center + width;
//		-|----------|------------|----------|-------
//		landa_i	  left         right    landa_i+1
		return (right - left) <= gapTol * right;
	}
	void shift(const double relShift,  const double alpha) {
		width += 2 * alpha * (fabs(center) + fabs(relShift));
		center += relShift;
	}
	bool isTooWide(const double relAcc) const {
		return width > relAcc * fabs(center);
	}
	bool isTooWide(const double relAcc, const double left, const double right) const {
		return isTooWide(relAcc) && left != center && right != center;
	}

	friend std::ostream& operator<<(std::ostream& os, const Interval & i) {
		os << i.center << " +/- " << i.width/2 << " " <<
			i.leftCount << ":" << i.rightCount << std::endl;
		return os;
	}
};


class GK  {// : public CompileTime {  //  GK matrix is a symmetric tridiagonal 
											//  square matrix of size NxN,
											//  having same entries on diagonal and 
											//  different offdiagonal 
public:
	enum Representation { 
		// in all representations indexing starts from 0
		Native = 0, // holds offdiagona original entries of GK matrix
						// T -> tridiagonal GK matrix ( with constant diagonal)
						// T(i,i) = data[0], for all i
						// T(i,i+1) = data[i+1]
					
		LDL = 1,    // T = LDL^T, i
						// L_i -> data[2 * i + 2]
						// D_i -> data[2 * i + 1]

		URU = 2,    // T = URU^T
						// shift ->data[0]
						// U_i -> data[2 * i + 2]
						// R_i -> data[2 * i + 1]
		
		NGN = 3, 	// T = NGN_k, for all k
						// L_i -> data[2 * i + 2]
						// D_i -> data[2 * i + 1]
						// U_i -> data[4 * n + 2 * i + 2]
						// R_i -> data[4 * n + 2 * i + 1]
						
		NGNPack = 4	// T = NGN_k, for single k
						// L_i -> data[2 * i + 2] , i < k
						// D_i -> data[2 * i + 1] , i < k
						// gamma_k -> data[2 * k + 1];
						// U_i -> data[2 * i + 2] , i > k
						// R_i -> data[2 * i + 1] , i > k
						// k is stored in repIndex
		// some comments
		// keep in mind that the following are equivalent
		// LDL <=> NGNPack with repIndex == 2 * n - 1
		// URU <=> NGNPack with repIndex == 0
	};

	// CONSTRUCTORS
	// create from BiDiagonal
	GK(const Matrix & B, const double shift_ = 0);
	GK(const Matrix & B, double * RS allocateData, const int space_, 
			const double shift_ = 0); 

	// Create from other GK	

	GK(const Representation rep, const GK & oth);
	GK(const Representation rep, const GK & oth, const double relShift_);

	GK(const Representation rep, const GK & oth,
		 double * RS allocateData, const int spase_, 
		const double relShift_ = 0);

	GK(const Representation rep, const int n_, 
		double * RS allocateData, const int spase_);

	GK(const Representation rep, const int n_);
	// empty constructor
	GK(){	selfAllocated = false; 	 sigAllocated = false; 	};

	void biDiagonalToNative(const double * RS d, const double * RS e, 
									const double shift);
	void biDiagonalToNative(const double * RS d, const double * RS e, 
									int * RS signature,  const double shift);

	virtual ~GK() { delocate(); };
	// static methods
	static int spaceNeeded(const int n, const Representation rep);
	static int spaceSignature(const int n){
		int signatureSize = 2 * n / sizeof(int);
		signatureSize += (2 * n ) % sizeof(int) > 0;
		return signatureSize;
	}

	// setters
	void setPointers();
	void setPivotMin(const double val) { pivotMin = val; };
	void setShift(const double shift_) { *this = shift_; }
	void shiftFrom(const GK & oth, const double relShift);
	void setMinGammaIndex(const int i) { minGammaIndex = i;};

	// getters
	int            getN              () const { return n; };
	double         getShift          () const { return *data; };
	Representation getRepresentation () const { return representation; };
	const double * getData           () const { return data; };
	const double * getNat            () const { return nat ; };
 	const double * getGamma          () const { return data + 8 * n; };
	const int    * getSig            () const { return sig ; };


	void getResidual    (const double * RS vec, const int incVec, 
								double * RS res      , const int incRes) const;
	double getPivotMin() const { return pivotMin; };
	const char * getRepresentationName() const;

	// SOLVERS
	void solveVec(const double * RS b, const int incB, double * RS x, const int incX) {		(this->*solve)(b, incB, x, incX);	}

	void getEigen(EIGEN_FUNC_PARAMS) const {
		(this->*getVec)(x, incX, xNorm2, residual, normUV);}

	
	void getMinGamma(int & k, double & gamma_k) const {
		k = minGammaIndex;
		if (isNGN()) 
			gamma_k = getGamma()[minGammaIndex];
		else gamma_k = -1;
	}
	int getMinGamma() const { return minGammaIndex; };
	int getRepIndex() const { 
		return repIndex >= 0 ? repIndex : chooseRepIndex(); };

	GK& operator+= (const double & d) { (this->*addShift)(d);   return *this;}
	GK& operator = (const double & d) { (this->*initialize)(d); return *this;}

	// when representation is computed, we want to know how accurate it is
	double checkRepresentation(int & k) const  {	return (this->*check)(k);	}
	double checkRepresentationNative (int & k) const;
	double checkRepresentationLDL    (int & k) const;
	double checkRepresentationURU    (int & k) const;
	double checkRepresentationNGN    (int & k) const;
	double checkRepresentationNGNPack(int & k) const;
	void checkRepresentationAllNGN(double * error) const;


	int countNative (const double relShift) const;
	int countLDL    (const double relShift) const;
	int countURU    (const double relShift) const;
	int countNGN    (const double relShift) const;
	int countNGNPack(const double relShift) const;
	int callCount   (const double relShift) const {
		return (this->*count)(relShift) - n;	}

	void refine(Interval & inter, const double relTol) const; 

	//comptes the local eigen values and stores them in intervals
//	void localEigen(Interval all[], const int clusterStart, 
//								const int clusterStop) const;

	double findEigenBorder(const int i, const double guess, 
								const double step) const;

	friend std::ostream& operator<<(std::ostream& os, const GK & m);
	void callPrint( std::ostream& os) const{	(this->*print)(os);	}
	// the following methods solve 
	// T x = b
	void solveVecNative (const double * RS b, const int incB, double * RS x, const int incX) const;
	void solveVecLDL    (const double * RS b, const int incB, double * RS x, const int incX) const;
	void solveVecURU    (const double * RS b, const int incB, double * RS x, const int incX) const;
	void solveVecNGN    (const double * RS b, const int incB, double * RS x, const int incX) const;
	void solveVecNGNPack (const double * RS b, const int incB, double * RS x, const int incX) const;

	// the following methods approximate eigen vector corresponding to 
	// eigen value 0
	// T x = 0 x
	// residual & norm are also calculated
	void getEigenLDL     (EIGEN_FUNC_PARAMS) const; 
	void getEigenURU     (EIGEN_FUNC_PARAMS) const; 
	void getEigenNGN     (EIGEN_FUNC_PARAMS) const; 
	void getEigenNGNPack (EIGEN_FUNC_PARAMS) const; 
	void getEigenNative  (EIGEN_FUNC_PARAMS) const; 
	void getEigenNative2 (EIGEN_FUNC_PARAMS) const; 
	void getEigenNative3 (EIGEN_FUNC_PARAMS) const; 


	// the following methods solve 
	// T x = e_k
	// where, e_k is the canonical vector having 1.0 at index k
	void solveCanonicalLDL   (const int k, EIGEN_FUNC_PARAMS) const;
	void solveCanonicalURU   (const int k, EIGEN_FUNC_PARAMS) const;
	void solveCanonicalNGN   (const int k, EIGEN_FUNC_PARAMS) const;
	void solveCanonicalNGN2  (const int k, EIGEN_FUNC_PARAMS) const;
	void solveCanonicalNGNPack(const int k, EIGEN_FUNC_PARAMS) const;
	void solveCanonicalNative(const int k, EIGEN_FUNC_PARAMS) const;

	int chooseRepIndex() const {
		int k;
		checkRepresentation(k); // find the best representation
		repIndex = k;
		return k;
	};
	void setRepIndex(const int choise) {repIndex = choise; };
	int chooseRepIndex(const double * RS s2, double & val);  //envelope vector

	void setNat(const double * val) { nat = val;};
	void setSig(const int    * val) { sig = val;};
	void setRepresentation(Representation newRep) {
		representation = newRep;
		setPointers();
	}
	
	
	bool isLDL    () const { return representation == LDL   ; };
	bool isURU    () const { return representation == URU   ; };
	bool isNGN    () const { return representation == NGN   ; };
	bool isNGNPack() const { return representation == NGNPack; };
	bool isNative () const { return representation == Native; };
	friend class GKManager;
private:
	void construct(const Representation rep, const int n_, double * RS allocateData, const int spase_);
   double getL    (const int i) const { return *(data + 2 + 2 * i); };
	double getD    (const int i) const { return *(data + 1 + 2 * i); };
	double getR    (const int i) const { return *(data + 4 * n + 1 + 2 * i); };
	double getU    (const int i) const { return *(data + 4 * n + 2 + 2 * i); };
	double getN    (const int i, const int k) const { 
		return i < k ? getL(i) : getU(i); };
	double getG    (const int i, const int k) const { 
		return i < k ? getD(i) : (i > k ? getR(i) : getGamma()[i]); };
	double getNative(const int i) const {
//		return (i & 1) == 0 ? e[(i >> 1) - 1] : d[i >> 1];
		assert(i >= 0 && i < 2 * n);
		return nat[i];
	}

	double pivot(const double d) const {
      return fabs(d) >= pivotMin ? d : -pivotMin;
	}
	double pivotLDL(const double d)  {
      return fabs(d) >= pivotMin ? d : -pivotMin;
	}
	void updateMinGamma(const double g, const int k, double & val) {	
		double tmp = fabs(g);
		minGammaIndex = tmp < val ? k   : minGammaIndex;
		val           = tmp < val ? tmp : val;
	}
	void copyFields(const GK & mat);
	double getShiftness(const double representationError, 
							  const double elementGrowth) const {
		double fact = elementGrowth < 8 * maxNat ? 0.5 : 1.0;
		return representationError + DBL_EPSILON * elementGrowth * fact;
	}


	
	
	

	void delocate();
	void nativeToLDL   (const GK & mat, const double shift);
	void nativeToURU   (const GK & mat, const double shift);
	void nativeToNGN   (const GK & mat, const double shift);
	void nativeToNGNPack  (const GK & mat, const double shift);
	void nativeToNative(const GK & mat, const double shift);

	void LDLToLDL      (const GK & mat, const double shift);
	void LDLToURU      (const GK & mat, const double shift);
	void LDLToNGN      (const GK & mat, const double shift);
	void LDLToNGNPack  (const GK & mat, const double shift);
	void LDLToNative   (const GK & mat, const double shift);

	void URUToLDL      (const GK & mat, const double shift);
	void URUToURU      (const GK & mat, const double shift);
	void URUToNGN      (const GK & mat, const double shift);
	void URUToNGNPack  (const GK & mat, const double shift);
	void URUToNative   (const GK & mat, const double shift);

	void NGNToLDL      (const GK & mat, const double shift);
	void NGNToURU      (const GK & mat, const double shift);
	void NGNToNGN      (const GK & mat, const double shift);
	void NGNToNGNPack  (const GK & mat, const double shift);
	void NGNToNative   (const GK & mat, const double shift);

	void NGNPackToNative   (const GK & mat, const double shift);
	void NGNPackToLDL      (const GK & mat, const double shift);
	void NGNPackToURU      (const GK & mat, const double shift);
	void NGNPackToNGN      (const GK & mat, const double shift);
	void NGNPackToNGNPack  (const GK & mat, const double shift);



	void initializeLDL   (const double shift);
	void initializeURU   (const double shift);
	void initializeNGN   (const double shift);
	void initializeNGNPack(const double shift);
	void initializeNative(const double shift);

	void printLDL    (std::ostream & os) const;
	void printURU    (std::ostream & os) const;
	void printNGN    (std::ostream & os) const;
	void printNGNPack(std::ostream & os) const;
	void printNative (std::ostream & os) const;

	void addShiftLDL    (const double & );
	void addShiftURU    (const double & );
	void addShiftNGN    (const double & );
	void addShiftNGNPack(const double & );
	void addShiftNative (const double & );

	factorize   						from[5];
	initializeF 						initialize;
	solveVecor  						solve; 
	getEigenVec 						getVec;
	printF      						print;
	addShiftF   						addShift;
	solveCanonicalVector 			canonicalSolve;
	checkRepresentationF				check;
	countF								count;



	double      * RS data; // 2n double represents  2 n - 1 degrees of freedom GK matrix + 1 shift :) 
	const double * RS nat; // pointer to native representation
	const int    * RS sig; // pointer to signature matrix
	double       pivotMin; // holds the pivot min threshold
	double			maxNat; // abs maximum offdiagonal element
	double 			spdiam; // estimate of diameter of the matrix
	int                 n; // 2 * n -> dimention of the GK matrix
	int             space; // external space provided to this class
	mutable int minGammaIndex; // index where |gamma_i| is min, if is NGN
	mutable int  repIndex; // holds good index to shift from
	Representation representation; // holds current representation
	bool    selfAllocated; // if the class allocated it's own mem
	bool     sigAllocated; // true, if this object allocated sig matrix
};
const GK::Representation  Native = GK::Native;
const GK::Representation  LDL    = GK::LDL   ;
const GK::Representation  URU    = GK::URU   ;
const GK::Representation  NGN    = GK::NGN   ;
const GK::Representation  NGNPack= GK::NGNPack;


class GKBisector {
public:
	// constructs the bisector
	GKBisector(const GK & gk);
	GKBisector(const GK & gk, double * RS storage);
	~GKBisector();
	void setGapTol(const double val) { gapTol = val; };
	void setRelAcc(const double val) { relAcc = val; };
	void setIntervalArray(Interval * arr) { all = arr; };
	// returns the number of singular values, less then lambda
	// complexity should be O(n)
	int count(const double lambda) const; 

	// refines a single eigenvalue to full accuracy
	// assumes that inter will contain only 1 eigen value
	void refine(Interval & inter) const; 
	// compute the intervals of singeltons and clusters
	//
	void refine(const int i) const; 
	int split(const int i) const ;

	void identifyClusters(Interval inter[] ) const;

	double getGapTol() const { return gapTol; };
	double getRelAcc() const { return relAcc; };
	Interval * getIntervalArray() { return all; };
	int getN() const { return n; };
	double getMax() const { return max; };
private:
	double pivot(const double d) const {
      return fabs(d) >= pivotMin ? d : -pivotMin;
	}
	void initialize(const GK & gk);	
	double * RS uru;	
	Interval * all;
	double gapTol; // gap tolerance
	double relAcc; // realite accuracy limit for intitial eigenvalues
	double relativeWidthLimit;
	double max;    
	double min;
	double pivotMin;
	int n; 
	bool selfAllocated;
};
